clc,clear
load('donors_father.mat');
donors_father = donors_father;
load('weight_donors.mat');
weight_donors = weight_donors;
folder_path = 'C:\Users\qx\Desktop\Training and testing datasets\Training Set';  % 替换为你的文件夹路径
[donors,acceptors] = extract_data(folder_path);
folder_path = 'C:\Users\qx\Desktop\Training and testing datasets\Training Set';
[donors_false,acceptors_false] = extract_false_data(folder_path);

j = 1;
for i = 1:length(donors)
    reg = regexp(donors(i,1),'[acgt]','match');
    if(length(reg{1}) == 18)
        donors_new(j,1) = donors(i);
        j = j+1;
    end
end
donors = donors_new;

[common,idx] = intersect(donors_false,donors);
donors_false(idx) = [];

bases = {'a', 'c', 'g', 't'};

scores = scores_calculate(donors,weight_donors,donors_father);
scores_sum = sum(scores,1);
scores_sum = sort(scores_sum,2,"descend");

scores_2 = scores_calculate(donors_false(3:5000),weight_donors,donors_father);
scores_sum_2 = sum(scores_2,1);
scores_sum_2 = sort(scores_sum_2,2,"descend");

j = 1;
for i = 12:0.1:30
    TP = scores_sum(find(scores_sum > i));
    FP = scores_sum_2(find(scores_sum_2 > i));
    TN = scores_sum_2(find(scores_sum_2 < i));
    FN = scores_sum(find(scores_sum < i));
    sn = length(TP)/(length(TP)+length(FN));
    sp = length(TN)/(length(TN)+length(FP));
    ac = length(TP)/(length(TP)+length(FP));
    precision = length(TP)/(length(TP)+length(FP));
    recall = length(TP)/(length(TP)+length(FN));
    f1(1,j) = i;
    f1(2,j) = 2*precision*recall/(precision+recall);
    AC(1,j) = i;
    AC(2,j) = ac;
    AUC(1,j) = sn;
    AUC(2,j) = sp;
    j = j+1;
end
%%% 阈值21
figure
plot(AUC(2,:),AUC(1,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('Sp-Sn Of donorclassifier')
ylabel('Sp')
xlabel('Sn')

figure
plot(AC(1,:),AC(2,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('Ac Of donorclassifier')
ylabel('AC')
xlabel('threshold')

figure
plot(1-AUC(2,:),AUC(1,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('ROC Of donorclassifier')
ylabel('Sn')
xlabel('FPR')
AUC_1 = trapz(1-AUC(2,:),AUC(1,:));

%%  acceptors
clc,clear;
load('acceptors_father.mat');
acceptors_father = acceptors_father;
load('weight_acceptors.mat');
weight_acceptors = weight_acceptors;
folder_path = 'C:\Users\qx\Desktop\Training and testing datasets\Training Set';  % 替换为你的文件夹路径
[acceptors,acceptors] = extract_data(folder_path);
folder_path = 'C:\Users\qx\Desktop\Training and testing datasets\Training Set';
[acceptors_false,acceptors_false] = extract_false_data(folder_path);

j = 1;
for i = 1:length(acceptors)
    reg = regexp(acceptors(i,1),'[acgt]','match');
    if(length(reg{1}) == 36)
        acceptors_new(j,1) = acceptors(i);
        j = j+1;
    end
end
acceptors = acceptors_new;

[common,idx] = intersect(acceptors_false,acceptors);
acceptors_false(idx) = [];

bases = {'a', 'c', 'g', 't'};

scores = scores_calculate_2(acceptors,weight_acceptors,acceptors_father);
scores_sum = sum(scores,1);
scores_sum = sort(scores_sum,2,"descend");

scores_2 = scores_calculate_2(acceptors_false,weight_acceptors,acceptors_father);
scores_sum_2 = sum(scores_2,1);
scores_sum_2 = sort(scores_sum_2,2,"descend");

j = 1;
for i = 22:0.1:50
    TP = scores_sum(find(scores_sum > i));
    FP = scores_sum_2(find(scores_sum_2 > i));
    TN = scores_sum_2(find(scores_sum_2 < i));
    FN = scores_sum(find(scores_sum < i));
    sn = length(TP)/(length(TP)+length(FN));
    sp = length(TN)/(length(TN)+length(FP));
    ac = length(TP)/(length(TP)+length(FP));
    precision = length(TP)/(length(TP)+length(FP));
    recall = length(TP)/(length(TP)+length(FN));
    f1(1,j) = i;
    f1(2,j) = 2*precision*recall/(precision+recall);
    AC(1,j) = i;
    AC(2,j) = ac;
    AUC(1,j) = sn;
    AUC(2,j) = sp;

    j = j+1;
end
%%% 阈值41.4
figure
plot(AUC(2,:),AUC(1,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('Sp-Sn Of acceptorclassifier')
ylabel('Sp')
xlabel('Sn')

figure
plot(AC(1,:),AC(2,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('Ac Of acceptorclassifier')
ylabel('AC')
xlabel('threshold')

figure
plot(1-AUC(2,:),AUC(1,:));
ax = gca;
set(gca,'TickDir','out');
ax.Box = 'off';
ax.YAxisLocation = 'left';  
ax.XAxisLocation = 'bottom'; 
title('ROC Of acceptorclassifier')
ylabel('Sn')
xlabel('FPR')
AUC_1 = trapz(1-AUC(2,:),AUC(1,:));